Accelerating pulsar timing data analysis

The analysis of pulsar timing data, especially in pulsar timing array (PTA) projects, has encountered practical difficulties: evaluating the likelihood and/or correlation-based statistics can become prohibitively computationally expensive for large datasets. In situations where a stochastic signal of interest has a power spectral density that dominates the noise in a limited bandwidth of the total frequency domain (e.g. the isotropic background of gravitational waves), a linear transformation exists that transforms the timing residuals to a basis in which virtually all the information about the stochastic signal of interest is contained in a small fraction of basis vectors. By only considering such a small subset of these "generalised residuals", the dimensionality of the data analysis problem is greatly reduced, which can cause a large speedup in the evaluation of the likelihood: the ABC-method (Acceleration By Compression). The compression fidelity, calculable with crude estimates of the signal and noise, can be used to determine how far a dataset can be compressed without significant loss of information. Both direct tests on the likelihood, and Bayesian analysis of mock data, show that the signal can be recovered as well as with an analysis of uncompressed data. In the analysis of IPTA Mock Data Challenge datasets, speedups of a factor of three orders of magnitude are demonstrated. For realistic PTA datasets the acceleration may become greater than six orders of magnitude due to the low signal to noise ratio

Gravitational wave detection and data analysis for pulsar timing arrays

Long-term precise timing of Galactic millisecond pulsars holds great promise for measuring long-period (months-to-years) astrophysical gravitational waves. In this work we develop a Bayesian data analysis method for projects called pulsar timing arrays; projects aimed to detect these gravitational waves by using millisecond pulsars as very precise extraterrestrial clocks. We also introduce the gravitational-wave memory effect, and we describe how to detect this effect with pulsar timing arrays with the Bayesian method. The Bayesian data analysis method is developed, extensively tested, and applied to real data. This has resulted in the most accurate upper-limit on the stochastic gravitational-wave background.UBL - phd migration 201

From bright binaries to bumpy backgrounds: Mapping realistic gravitational wave skies with pulsar-timing arrays

Within the next several years, pulsar-timing array programs will likely usher in the next era of gravitational-wave astronomy through the detection of a stochastic background of nanohertz-frequency gravitational waves, originating from a cosmological population of inspiraling supermassive binary black holes. While the source positions will likely be isotropic to a good approximation, the gravitational-wave angular power distribution will be anisotropic, with the most massive and/or nearby binaries producing signals that may resound above the background. We study such a realistic angular power distribution, developing fast and accurate sky-mapping strategies to localize pixels and extended regions of excess power while simultaneously modeling the background signal from the less massive and more distant ensemble. We find that power anisotropy will be challenging to discriminate from isotropy for realistic gravitational-wave skies, requiring SNR >10 in order to favor anisotropy with 10:1 posterior odds in our case study. Amongst our techniques, modeling the population signal with multiple point sources in addition to an isotropic background provides the most physically motivated and easily interpreted maps, while spherical-harmonic modeling of the square-root power distribution, P(Ω)^(1/2), performs best in discriminating from overall isotropy. Our techniques are modular and easily incorporated into existing pulsar-timing array analysis pipelines

The Need For Speed: Rapid Refitting Techniques for Bayesian Spectral Characterization of the Gravitational Wave Background Using PTAs

Current pulsar timing array (PTA) techniques for characterizing the spectrum of a nanohertz-frequency stochastic gravitational-wave background (SGWB) begin at the stage of timing data. This can be slow and memory intensive with computational scaling that will worsen PTA analysis times as more pulsars and observations are added. Given recent evidence for a common-spectrum process in PTA data sets and the need to understand present and future PTA capabilities to characterize the SGWB through large-scale simulations, we have developed efficient and rapid approaches that operate on intermediate SGWB analysis products. These methods refit SGWB spectral models to previously-computed Bayesian posterior estimations of the timing power spectra. We test our new methods on simulated PTA data sets and the NANOGrav $12.5$-year data set, where in the latter our refit posterior achieves a Hellinger distance from the current full production-level pipeline that is $\lesssim 0.1$. Our methods are $\sim10^2$--$10^4$ times faster than the production-level likelihood and scale sub-linearly as a PTA is expanded with new pulsars or observations. Our methods also demonstrate that SGWB spectral characterization in PTA data sets is driven by the longest-timed pulsars with the best-measured power spectral densities which is not necessarily the case for SGWB detection that is predicated on correlating many pulsars. Indeed, the common-process spectral properties found in the NANOGrav $12.5$-year data set are given by analyzing only the $\sim10$ longest-timed pulsars out of the full $45$ pulsar array, and we find that the ``shallowing'' of the common-process power-law model occurs when gravitational-wave frequencies higher than $\sim 50$~nanohertz are included. The implementation of our methods is openly available as a software suite to allow fast and flexible PTA SGWB spectral characterization and model selection.Comment: 19 pages, 12 figures. Submitting to Physical Review

Efficient Gravitational Wave Searches with Pulsar Timing Arrays using Hamiltonian Monte Carlo

Pulsar timing arrays (PTAs) detect low-frequency gravitational waves (GWs) bylooking for correlated deviations in pulse arrival times. Current Bayesiansearches use Markov Chain Monte Carlo (MCMC) methods, which struggle to samplethe large number of parameters needed to model the PTA and GW signals. As thedata span and number of pulsars increase, this problem will only worsen. Analternative Monte Carlo sampling method, Hamiltonian Monte Carlo (HMC),utilizes Hamiltonian dynamics to produce sample proposals informed byfirst-order gradients of the model likelihood. This in turn allows it toconverge faster to high dimensional distributions. We implement HMC as analternative sampling method in our search for an isotropic stochastic GWbackground, and show that this method produces equivalent statistical resultsto similar analyses run with standard MCMC techniques, while requiring 100-200times fewer samples. We show that the speed of HMC sample generation scales as$\mathcal{O}(N_\mathrm{psr}^{5/4})$ where $N_\mathrm{psr}$ is the number ofpulsars, compared to $\mathcal{O}(N_\mathrm{psr}^2)$ for MCMC methods. Thesefactors offset the increased time required to generate a sample using HMC,demonstrating the value of adopting HMC techniques for PTAs.<br

A data analysis library for gravitational wave detection

One of the main goals of Pulsar Timing Arrays (PTAs) is the direct detection of gravitational waves (GWs). A first detection will be a major leap for astronomy and substantial effort is currently going into timing as many pulsars as possible, with the highest possible accuracy. As part of the individual PTA projects, several groups are developing data analysis methods for the final stage of a gravitational-waves search pipeline: the analysis of the timing residuals. Here we report the progress of on-going work to develop, within a Bayesian framework, a comprehensive and user friendly analysis library to search for gravitational waves in PTA data

Analysis of the first IPTA Mock Data Challenge by the EPTA timing data analysis working group

This is a summary of the methods we used to analyse the first IPTA Mock Data Challenge (MDC), and the obtained results. We have used a Bayesian analysis in the time domain, accelerated using the recently developed ABC-method which consists of a form of lossy linear data compression. The TOAs were first processed with Tempo2, where the design matrix was extracted for use in a subsequent Bayesian analysis. We used different noise models to analyse the datasets: no red noise, red noise the same for all pulsars, and individual red noise per pulsar. We sampled from the likelihood with four different samplers: "emcee", "t-walk", "Metropolis-Hastings", and "pyMultiNest". All but emcee agreed on the final result, with emcee failing due to artefacts of the high-dimensionality of the problem. An interesting issue we ran into was that the prior of all the 36 (red) noise amplitudes strongly affects the results. A flat prior in the noise amplitude biases the inferred GWB amplitude, whereas a flat prior in log-amplitude seems to work well. This issue is only apparent when using a noise model with individually modelled red noise for all pulsars. Our results for the blind challenges are in good agreement with the injected values. For the GWB amplitudes we found h_c = 1.03 +/- 0.11 [10^{-14}], h_c = 5.70 +/- 0.35 [10^{-14}], and h_c = 6.91 +/- 1.72 [10^{-15}], and for the GWB spectral index we found gamma = 4.28 +/- 0.20, gamma = 4.35 +/- 0.09, and gamma = 3.75 +/- 0.40. We note that for closed challenge 3 there was quite some covariance between the signal and the red noise: if we constrain the GWB spectral index to the usual choice of gamma = 13/3, we obtain the estimates: h_c = 10.0 +/- 0.64 [10^{-15}], h_c = 56.3 +/- 2.42 [10^{-15}], and h_c = 4.83 +/- 0.50 [10^{-15}], with one-sided 2 sigma upper-limits of: h_c <= 10.98 [10^{-15}], h_c <= 60.29 [10^{-15}], and h_c <= 5.65 [10^{-15}]

Model-based asymptotically optimal dispersion measure correction for pulsar timing

In order to reach the sensitivity required to detect gravitational waves, pulsar timing array experiments need to mitigate as much noise as possible in timing data. A dominant amount of noise is likely due to variations in the dispersion measure. To correct for such variations, we develop a statistical method inspired by the maximum likelihood estimator and optimal filtering. Our method consists of two major steps. First, the spectral index and amplitude of dispersion measure variations are measured via a time-domain spectral analysis. Second, the linear optimal filter is constructed based on the model parameters found in the first step, and is used to extract the dispersion measure variation waveforms. Compared to current existing methods, this method has better time resolution for the study of short timescale dispersion variations, and generally produces smaller errors in waveform estimations. This method can process irregularly sampled data without any interpolation because of its time-domain nature. Furthermore, it offers the possibility to interpolate or extrapolate the waveform estimation to regions where no data is available. Examples using simulated data sets are included for demonstration.Comment: 15 pages, 15 figures, submitted 15th Sept. 2013, accepted 2nd April 2014 by MNRAS. MNRAS, 201

European Pulsar Timing Array limits on continuous gravitational waves from individual supermassive black hole binaries

We have searched for continuous gravitational wave (CGW) signals produced by individually resolvable, circular supermassive black hole binaries (SMBHBs) in the latest European Pulsar Timing Array (EPTA) data set, which consists of ultraprecise timing data on 41-ms pulsars. We develop frequentist and Bayesian detection algorithms to search both for monochromatic and frequency-evolving systems. None of the adopted algorithms show evidence for the presence of such a CGW signal, indicating that the data are best described by pulsar and radiometer noise only. Depending on the adopted detection algorithm, the 95 per cent upper limit on the sky-averaged strain amplitude lies in the range 6 × 10^(−15) 10^9M_⊙ out to a distance of about 25 Mpc, and with M_c>10^(10)M_⊙ out to a distance of about 1Gpc (z ≈ 0.2). We show that state-of-the-art SMBHB population models predict <1 per cent probability of detecting a CGW with the current EPTA data set, consistent with the reported non-detection. We stress, however, that PTA limits on individual CGW have improved by almost an order of magnitude in the last five years. The continuing advances in pulsar timing data acquisition and analysis techniques will allow for strong astrophysical constraints on the population of nearby SMBHBs in the coming years

Are we there yet? Time to detection of nanohertz gravitational waves based on pulsar-timing array limits

Decade-long timing observations of arrays of millisecond pulsars have placed highly constraining upper limits on the amplitude of the nanohertz gravitational-wave stochastic signal from the mergers of supermassive black hole binaries (~10^(−15) strain at f = 1 yr^(−1)). These limits suggest that binary merger rates have been overestimated, or that environmental influences from nuclear gas or stars accelerate orbital decay, reducing the gravitational-wave signal at the lowest, most sensitive frequencies. This prompts the question whether nanohertz gravitational waves (GWs) are likely to be detected in the near future. In this Letter, we answer this question quantitatively using simple statistical estimates, deriving the range of true signal amplitudes that are compatible with current upper limits, and computing expected detection probabilities as a function of observation time. We conclude that small arrays consisting of the pulsars with the least timing noise, which yield the tightest upper limits, have discouraging prospects of making a detection in the next two decades. By contrast, we find large arrays are crucial to detection because the quadrupolar spatial correlations induced by GWs can be well sampled by many pulsar pairs. Indeed, timing programs that monitor a large and expanding set of pulsars have an ~80% probability of detecting GWs within the next 10 years, under assumptions on merger rates and environmental influences ranging from optimistic to conservative. Even in the extreme case where 90% of binaries stall before merger and environmental coupling effects diminish low-frequency gravitational-wave power, detection is delayed by at most a few years